library(ggordiplots)
library(microViz)
library("vegan")
library(phyloseqCompanion)
library(ggplot2)
library(Rmisc)
library(ggsignif)
library(ggpubr)
library(dplyr)
library("gridExtra")
library(microbiomeMarker)
library(mblm)
library(tidyr)
require(MASS)
require(boot)
require(pscl)

# set directory
setwd("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing/Code")

# set scientific notation off
options(scipen=999)

# load exploration type and guild data
load("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing/Code/RA otu by guild and exploration type.RData")

soil.data <- read.csv("soil.data.final.csv")

# load carcass data
sample_data(physeq)$decomposing.carcass <- soil.data$decomposing.carcass[1:96]
sample_data(physeq)$organic.soil.C.N <- soil.data$organic.soil.C.N[1:96]
sample_data(physeq)$organic.soil.GWC <- soil.data$organic.soil.GWC[1:96]
sample_data(physeq)$carcass.deposition <- soil.data$carcass.deposition[1:96]
sample_data(physeq)$ph <- soil.data$ph[1:96]
sample_data(physeq)$white.spruce <- soil.data$white.spruce[1:96]
sample_data(physeq)$paper.birch <- soil.data$paper.birch[1:96]
sample_data(physeq)$dwarf.birch <- soil.data$dwarf.birch[1:96]
sample_data(physeq)$moss <- soil.data$moss[1:96]
sample_data(physeq)$vaccinium <- soil.data$vaccinium[1:96]
sample_data(physeq)$fern <- soil.data$fern[1:96]
sample_data(physeq)$kinnikinnick <- soil.data$kinnikinnick[1:96]
sample_data(physeq)$heather <- soil.data$heather[1:96]
sample_data(physeq)$eriophorum <- soil.data$eriophorum[1:96]
sample_data(physeq)$fireweed <- soil.data$fireweed[1:96]
sample_data(physeq)$horsetail <- soil.data$horsetail[1:96]
sample_data(physeq)$alder <- soil.data$alder[1:96]
sample_data(physeq)$hogsweed <- soil.data$hogsweed[1:96]
sample_data(physeq)$willow <- soil.data$willow[1:96]
sample_data(physeq)$redberry <- soil.data$redberry[1:96]
sample_data(physeq)$azalea <- soil.data$azalea[1:96]
sample_data(physeq)$slope <- soil.data$slope[1:96]

# filter phyloseq object by Hansen stream
physeq_close <- subset_samples(physeq, Distance < 11)
physeq_final <- subset_samples(physeq_close, Stream=="Hansen")
physeq_final <- subset_samples(physeq_close, Stream=="Yako")
physeq_final <- subset_samples(physeq_close, Stream=="Happy")


# get numbers for alpha diversity metrics at Hansen
alpha.diversity <- estimate_richness(physeq, measures=c("Observed","Shannon", "Simpson"))
head(alpha.diversity)

# combine phyloseq object with alpha diversity metrics 
data <- cbind(sample_data(physeq_final), alpha.diversity)

# if not normally distributed, do unpaired two-samples Wilcoxon test
yes.observed <- data$Observed[which(data$decomposing.carcass=="yes")]
no.observed <- data$Observed[which(data$decomposing.carcass=="no")]
t.test(yes.observed, no.observed)
wilcox.test(yes.observed, no.observed, alternative = "two.sided", exact=F)

# at Happy, higher species richness with no carcass

# if not normally distributed, do unpaired two-samples Wilcoxon test
yes.diversity <- data$Shannon[which(data$decomposing.carcass=="yes")]
no.diversity <- data$Shannon[which(data$decomposing.carcass=="no")]
t.test(yes.diversity, no.diversity)
wilcox.test(yes.diversity, no.diversity, alternative = "two.sided", exact=F)

# make sure all four locations for exploration.guild.RA have "yes" for decomposing carcass (only for Hansen)
exploration.guild.RA$decomposing.carcass[which(exploration.guild.RA$Stream=="Hansen" & exploration.guild.RA$Side=="SR" & exploration.guild.RA$Distance==3 & 
                                                 exploration.guild.RA$Transect=="S4")]

exploration.guild.RA$decomposing.carcass[which(exploration.guild.RA$Stream=="Hansen" & exploration.guild.RA$Side=="SR" & exploration.guild.RA$Distance==1 & 
                                                 exploration.guild.RA$Transect=="S4")] <- "yes"

exploration.guild.RA$decomposing.carcass[which(exploration.guild.RA$Stream=="Hansen" & exploration.guild.RA$Side=="SR" & exploration.guild.RA$Distance==6 & 
                                                 exploration.guild.RA$Transect=="S4")]

exploration.guild.RA$decomposing.carcass[which(exploration.guild.RA$Stream=="Hansen" & exploration.guild.RA$Side=="SL" & exploration.guild.RA$Distance==10 & 
                                                 exploration.guild.RA$Transect=="S4")]

################################################################
# test richness of fungal guilds (sapros or symbios) with decomposing carcass
################################################################

exploration.guild.RA <- exploration.guild.RA[which(exploration.guild.RA$Stream=="Happy"),]
exploration.guild.RA <- exploration.guild.RA[!(exploration.guild.RA$Transect == "Bog"),]
exploration.guild.RA <- exploration.guild.RA[!(exploration.guild.RA$Stream == "Yako" & exploration.guild.RA$Transect =="Y3" & exploration.guild.RA$Distance == 10),]
exploration.guild.RA <- exploration.guild.RA[!(exploration.guild.RA$Stream == "Yako" & exploration.guild.RA$Transect =="Y2" & exploration.guild.RA$Distance == 10),]
exploration.guild.RA <- exploration.guild.RA[!(exploration.guild.RA$Stream == "Happy" & exploration.guild.RA$Transect =="H2" & exploration.guild.RA$Distance == 1),]


yes <- exploration.guild.RA$saprotroph.richness[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="yes")]
no <- exploration.guild.RA$saprotroph.richness[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="no")]
t.test(yes,no)
wilcox.test(yes,no, exact=F)

# saprotrophic richness higher with decomposing carcasses (Hansen and overall)

yes <- exploration.guild.RA$symbiotroph.richness[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="yes")]
no <- exploration.guild.RA$symbiotroph.richness[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="no")]
t.test(yes, no)
wilcox.test(yes, no, exact=F)

# at Happy, symbiotrophic richness was ALMOST higher with no carcass

################################################################
# test richness of exploration types with decomposing carcass
################################################################

yes <- exploration.guild.RA$short.richness[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="yes")]
no <- exploration.guild.RA$short.richness[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="no")]
t.test(yes, no)
wilcox.test(yes, no, exact=F)

# with new exploration.guild.RA data, short.richness higher with carcass!
# with new exploration guild RA data, short richness higher wiht NO carcass at Happy!

yes <- exploration.guild.RA$long.richness[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="yes")]
no <- exploration.guild.RA$long.richness[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="no")]
t.test(yes, no)
wilcox.test(yes, no, exact=F)

# no effect

################################################################
# test diversity of guilds with decomposing carcass
################################################################

yes <- exploration.guild.RA$saprotroph.diversity[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="yes")]
no <- exploration.guild.RA$saprotroph.diversity[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="no")]
t.test(yes, no)
wilcox.test(yes, no, exact=F)

# saprotrophic diversity was higher at decomposing carcasses (overall and Hansen)

yes <- exploration.guild.RA$symbiotroph.diversity[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="yes")]
no <- exploration.guild.RA$symbiotroph.diversity[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="no")]
t.test(yes, no)
wilcox.test(yes, no, exact=F)

# no effect

################################################################
# test diversity of exploration types with decomposing carcass
################################################################

yes <- exploration.guild.RA$short.diversity[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="yes")]
no <- exploration.guild.RA$short.diversity[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="no")]
t.test(yes, no)
wilcox.test(yes, no, exact=F)


yes <- exploration.guild.RA$long.diversity[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="yes")]
no <- exploration.guild.RA$long.diversity[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="no")]
t.test(yes, no)
wilcox.test(yes, no, exact=F)

# long diversity was higher with no decomposing carcass at Hansen
# long diversity had no effect with decomposing carcass with all data

#########################################
# test RA of fungal guilds with decomposing carcass
#########################################

left <- exploration.guild.RA$saprotroph[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="yes")]
right <- exploration.guild.RA$saprotroph[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="no")]
t.test(left, right)
wilcox.test(left, right, exact=F)

left <- exploration.guild.RA$symbiotroph[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="yes")]
right <- exploration.guild.RA$symbiotroph[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="no")]
t.test(left, right)
wilcox.test(left, right, exact=F)

# symbiotroph abundance is higher with carcass at Hansen
# symbiotroph abundance did not vary with carcass for all data

################################################################
# test RA of exploration types with decomposing carcass
################################################################

hansen.left.close <- exploration.guild.RA$medium.distance[which(exploration.guild.RA$decomposing.carcass=="yes")]
hansen.right.close <- exploration.guild.RA$medium.distance[which(exploration.guild.RA$decomposing.carcass=="no")]
t.test(hansen.left.close, hansen.right.close)
wilcox.test(hansen.left.close, hansen.right.close, exact=F)

# short distance relative abundance increased with carcass at Hansen
# short distsance did not vary with carcass with all data

hansen.left.close <- exploration.guild.RA$medium.distance[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="yes")]
hansen.right.close <- exploration.guild.RA$medium.distance[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="no")]
t.test(hansen.left.close, hansen.right.close)
wilcox.test(hansen.left.close, hansen.right.close, exact=F)

# at Happy, medium distance was higher with no carcass

hansen.left.close <- exploration.guild.RA$long.distance[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="yes" & exploration.guild.RA$Stream=="Hansen")]
hansen.right.close <- exploration.guild.RA$long.distance[which(exploration.guild.RA$Distance < 11 & exploration.guild.RA$decomposing.carcass=="no" & exploration.guild.RA$Stream=="Hansen")]
t.test(hansen.left.close, hansen.right.close)
wilcox.test(hansen.left.close, hansen.right.close, exact=F)

# no effect

################################################################
# BETA DIVERSITY with decomposing carcass
################################################################

physeq.ord <- ordinate(physeq_close, "PCoA", "bray")
physeq.ord <- ordinate(physeq_close, "NMDS", "bray")

# ANOSIM
transect_group = get_variable(physeq_close, "decomposing.carcass")
transect_ano = anosim(distance(physeq_close, "bray"), transect_group, permutations = 9999)
transect_ano$signif
transect_ano$statistic

# not significant difference in communities

##################################################################################################
# Wilcoxon test to compare the abundance of different genera with presence of decomposing carcass
####################################################################################################

# use relative abundance
physeq.prop <- transform_sample_counts(physeq, function(x) x/sum(x))

# subset samples
physeq_nobog <- subset_samples(physeq.prop, Transect != "Bog")

# subset by Hansen and agglomerate by genus
#physeq_hansen <- subset_samples(physeq.prop, Stream == "Hansen")
#physeq_close <- subset_samples(physeq.prop, Distance < 11)

# agglomerate by genus
genus.glom <- tax_glom(physeq_nobog, taxrank = 'Genus')

# get genus abundance matrix from phyloseq object
dat_genus <- psmelt(genus.glom)
dat_genus$Genus <- as.character(dat_genus$Genus)
dat_genus$Genus <- dat_genus$Genus %>% replace_na('Unclassified')
genus.abund <- subset(dat_genus, select = c("Sample", "Abundance", "Stream", "Transect", "Side", "Distance", "organic.soil.C.N", "organic.soil.GWC",
                                            "carcass.deposition", "ph", "decomposing.carcass", "Genus", "white.spruce", "paper.birch",
                                            "dwarf.birch", "moss", "vaccinium", "fern", "kinnikinnick", "heather", "eriophorum",
                                            "fireweed", "horsetail", "alder", "hogsweed", "willow", "redberry", "azalea","slope"))

# create unique number for each genera
genus.abund <- genus.abund %>%
  group_by(Genus) %>%
  mutate(ID = cur_group_id())

genus.abund <- as.data.frame(genus.abund)

# create vector for each genera and p value
pvalue <- c(1:308)
value <- c(1:308)

# or can just run wilcox test
for (i in 1:308){
  genus.final <- genus.abund[which(genus.abund$ID==i),]
  yes <- genus.final$Abundance[which(genus.final$decomposing.carcass=="yes")]
  no <- genus.final$Abundance[which(genus.final$decomposing.carcass=="no")]
  w <- wilcox.test(yes, no, exact=F)
  ttest <- t.test(yes, no)
  if (ttest$estimate[1] > ttest$estimate[2]){
    value[i] <- 1
  } else {value[i] <- 2}
  pvalue[i] <- w$p.value
}

# create dataframe with genus and p value combined
genus.p <- data.frame(levels(factor(genus.abund$Genus)), pvalue, value)
genus.p <- genus.p[order(pvalue),]
genus.sig <- genus.p[which(genus.p$pvalue < 0.051),]
colnames(genus.sig) <- c("Genus","P value")
View(genus.sig)

# assign lifestyle to each genera from above
# upload funguild database
funguild <- read.csv("Fungal Functional Traits database.csv", header=TRUE)
genus.sig$lifestyle <- c(1:85)
genus.sig$type <- c(1:85)
for (i in 1:85){
  genera <- genus.sig$Genus[i]
  genus.sig$lifestyle[i] <- funguild$primary_lifestyle[which(funguild$GENUS==genera)]
  genus.sig$type[i] <- funguild$Ectomycorrhiza_exploration_type_template[which(funguild$GENUS==genera)]
}

# export table of these genera
write.csv(genus.sig, "Genera with higher RA at decomposing carcass.csv")

# can run glm or lm for each genera
p.mod.carcass <- rep(0, 308)
p.mod.deposition <- rep(0, 308)
value.mod.deposition <- rep(0, 308)
value.mod.carcass <- rep(0, 308)
mean.abundance <- rep(0, 308)

for (i in 1:308){
  tryCatch({
  genus.final <- genus.abund[which(genus.abund$ID==i),]
  # run linear model with all variables and select best model using AIC before withdrawing
  # significance of individual predictors
  all.parms <- lm(Abundance ~ organic.soil.C.N + organic.soil.GWC + decomposing.carcass + carcass.deposition + Distance + white.spruce + paper.birch + dwarf.birch + moss + vaccinium + fern + kinnikinnick + 
              heather + eriophorum + fireweed + horsetail + alder + hogsweed + willow + redberry + azalea + slope, data=genus.final)
  step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
  temp <- summary(step.model)
  bleh <- temp$coefficients
  bleh <- as.data.frame(bleh)
  bleh <- tibble::rownames_to_column(bleh, "Variable")
  mean.abundance[i] <- mean(genus.final$Abundance)
  p.mod.deposition[i] <- bleh$`Pr(>|t|)`[which(bleh$Variable=="carcass.depositionyes")]
  value.mod.deposition[i] <- bleh$`Estimate`[which(bleh$Variable=="carcass.depositionyes")]
  p.mod.carcass[i] <- bleh$`Pr(>|t|)`[which(bleh$Variable=="decomposing.carcassyes")]
  value.mod.carcass[i] <- bleh$`Estimate`[which(bleh$Variable=="decomposing.carcassyes")]
  # just run linear model
  #mod <- glm.nb(Abundance ~ decomposing.carcass, data=genus.final)
  #mod <- lm(Abundance ~ organic.soil.C.N + organic.soil.GWC + decomposing.carcass + carcass.deposition + Distance + white.spruce + paper.birch + dwarf.birch + moss + vaccinium + fern + kinnikinnick + 
  #                  heather + eriophorum + fireweed + horsetail + alder + hogsweed + willow + redberry + azalea + slope, data=genus.final)
  #temp <- summary(mod)
  #p.mod.carcass[i] <- temp$coefficients[4,4]
  #p.mod.deposition[i] <- temp$coefficients[5,4]]) 
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
}

# CARCASS DEPOSITION # 
# create dataframe with genus and p value combined for carcass deposition
genus.p <- data.frame(levels(factor(genus.abund$Genus)), p.mod.deposition, value.mod.deposition, mean.abundance)
genus.p <- genus.p[order(p.mod.deposition),]
genus.sig <- genus.p[which(genus.p$p.mod.deposition > 0),]
colnames(genus.sig) <- c("Genus","P value", "Estimate")
View(genus.sig)

# assign lifestyle to each genera from above
# upload funguild database
funguild <- read.csv("Fungal Functional Traits database.csv", header=TRUE)
genus.sig$lifestyle <- c(1:53)
genus.sig$type <- c(1:53)
for (i in 1:53){
  genera <- genus.sig$Genus[i]
  genus.sig$lifestyle[i] <- funguild$primary_lifestyle[which(funguild$GENUS==genera)]
  genus.sig$type[i] <- funguild$Ectomycorrhiza_exploration_type_template[which(funguild$GENUS==genera)]
}

# export table of these genera
write.csv(genus.sig, "Genera with higher RA at carcass deposition.csv")

# DECOMPOSING CARCASS # 
# create dataframe with genus and p value combined for decomposing carcass
genus.p <- data.frame(levels(factor(genus.abund$Genus)), p.mod.carcass, value.mod.carcass, mean.abundance)
genus.p <- genus.p[order(p.mod.carcass),]
genus.sig <- genus.p[which(genus.p$p.mod.carcass > 0),]
colnames(genus.sig) <- c("Genus","P value", "Estimate")
View(genus.sig)

# assign lifestyle to each genera from above
# upload funguild database
funguild <- read.csv("Fungal Functional Traits database.csv", header=TRUE)
genus.sig$lifestyle <- c(1:8)
genus.sig$type <- c(1:8)
for (i in 1:8){
  genera <- genus.sig$Genus[i]
  genus.sig$lifestyle[i] <- funguild$primary_lifestyle[which(funguild$GENUS==genera)]
  genus.sig$type[i] <- funguild$Ectomycorrhiza_exploration_type_template[which(funguild$GENUS==genera)]
}

# export table of these genera
write.csv(genus.sig, "Genera with higher RA at decomposing carcass.csv")

# examine genera manually
genus.final <- genus.abund[which(genus.abund$Genus=="Xerocomus"),]
yes <- genus.final$Abundance[which(genus.final$decomposing.carcass=="yes")]
no <- genus.final$Abundance[which(genus.final$decomposing.carcass=="no")]

hansen <- genus.final[which(genus.final$Stream=="Hansen"),]
all.parms <- lm(Abundance ~ organic.soil.C.N + decomposing.carcass + carcass.deposition + Distance + white.spruce + paper.birch + dwarf.birch + moss + vaccinium + fern + kinnikinnick + 
                   heather + eriophorum + fireweed + horsetail + alder + hogsweed + willow + redberry + azalea + slope, data=genus.final)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

t.test(yes, no)
wilcox.test(yes, no, exact=F)

# export means by stream for each genus to make bar graph in excel
final <- genus.abund[genus.abund$Genus %in% genus.sig$Genus,]
tgc <- summarySE(final, measurevar="Abundance", groupvars=c("Genus","decomposing.carcass"))
tgc_order <- tgc[order(-tgc$Abundance),]
tgc_order
#write.csv(tgc, "means of sig genus abundance by Hansen side.csv")

#####################
# genera only found at decomposing carcass # 
#########################

# GENERA ONLY FOUND AT DECOMPOSING CARCASS SITES # 
sum.abund <- c(1:308)
for (i in 1:308){
    genus.final <- unique(genus.abund$Genus)[i]
    no.carcass <- genus.abund[which(genus.abund$decomposing.carcass=="no" & genus.abund$Genus==genus.final),]
    sum.abund[i] <- sum(no.carcass$Abundance)
}

test <- data.frame(unique(genus.abund$Genus), sum.abund)
colnames(test) <- c("Genus", "Sum")
View(test)

# assign functional guilds to each of the final genera
funguild <- read.csv("Fungal Functional Traits database.csv", header=TRUE)
test$lifestyle <- c(1:308)
test$type <- c(1:308)
for (i in 1:308){
  genera <- test$Genus[i]
  test$lifestyle[i] <- funguild$primary_lifestyle[which(funguild$GENUS==genera)]
  test$type[i] <- funguild$Ectomycorrhiza_exploration_type_template[which(funguild$GENUS==genera)]
}

test.final <- test[which(test$Sum == 0),]
View(test.final)

write.csv(test, "test.csv")

# SPECIES ONLY FOUND AT DECOMPOSING CARCASS SITES # 
sum.abund <- c(1:531)
for (i in 1:531){
  species.final <- unique(species.abund$Species)[i]
  no.carcass <- species.abund[which(species.abund$decomposing.carcass=="no" & species.abund$Species==species.final),]
  sum.abund[i] <- sum(no.carcass$Abundance)
}

test <- data.frame(unique(species.abund$Species), sum.abund)
colnames(test) <- c("Species", "Sum")
View(test)

# assign functional guilds to each of the final genera
funguild <- read.csv("Fungal Functional Traits database.csv", header=TRUE)
test$lifestyle <- c(1:531)
test$type <- c(1:531)
for (i in 1:531){
  genera <- test$Species[i]
  test$lifestyle[i] <- funguild$primary_lifestyle[which(funguild$GENUS==genera)]
  test$type[i] <- funguild$Ectomycorrhiza_exploration_type_template[which(funguild$GENUS==genera)]
}

test.final <- test$Species[which(test$Sum == 0)]

genera.final <- c(1:69)
for (i in 1:69){
  species.final <- test.final[i]
  genera.final[i] <- unique(all.otu.abundance.final$Genus[which(all.otu.abundance.final$Species==species.final)])
}

genera.final <- as.data.frame(genera.final)
genera.final$Species <- test.final
colnames(genera.final) <- c("Genus", "Species")

# assign functional guilds to each of the final genera
funguild <- read.csv("Fungal Functional Traits database.csv", header=TRUE)
genera.final$lifestyle <- c(1:69)
genera.final$type <- c(1:69)
for (i in 1:69){
  genera <- genera.final$Genus[i]
  genera.final$lifestyle[i] <- funguild$primary_lifestyle[which(funguild$GENUS==genera)]
  genera.final$type[i] <- funguild$Ectomycorrhiza_exploration_type_template[which(funguild$GENUS==genera)]
}

write.csv(genera.final, "species only at decomposing carcass.csv")

# take all genera that are at decomposing carcass
genera.carcass <- unique(genus.abund$Genus[which(genus.abund$decomposing.carcass=="yes" & genus.abund$Abundance > 0)])

genera.nocarcass <- unique(genus.abund$Genus[which(genus.abund$decomposing.carcass=="no" & genus.abund$Abundance == 0)])

intersect <- intersect(genera.carcass, genera.nocarcass)

genera.only.carcass <- setdiff(genera.carcass, genera.nocarcass)


species.carcass <- unique(species.abund$Species[which(species.abund$decomposing.carcass=="yes" & species.abund$Abundance > 0)])

species.nocarcass <- unique(species.abund$Species[which(species.abund$decomposing.carcass=="no" & species.abund$Abundance > 0)])

intersect <- intersect(species.carcass, species.nocarcass)

species.only.carcass <- setdiff(species.carcass, species.nocarcass)

genera.final <- c(1:66)
for (i in 1:66){
  species.final <- species.only.carcass[i]
  genera.final[i] <- unique(all.otu.abundance.final$Genus[which(all.otu.abundance.final$Species==species.final)])
}


####################################
# make figures # 
####################################

# just at Hansen
temp <- exploration.guild.RA[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 11),]

# SAPROTROPH RICHNESS # 

tgc <- summarySE(exploration.guild.RA, measurevar="saprotroph.richness", groupvars=c("decomposing.carcass"))
tgc

saprotroph.richness <- ggplot(tgc, aes(x=reorder(decomposing.carcass, -saprotroph.richness), y= saprotroph.richness, fill=decomposing.carcass)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() +
  theme_classic() + ylab("Species richness") + geom_errorbar(aes(ymin=saprotroph.richness-se, ymax=saprotroph.richness+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.y = element_text(size = 18), axis.text.y = element_text(size = 18), axis.text.x = element_text(size = 18), axis.title.x = element_text(size = 18),
        legend.position="none", plot.title = element_text(size=20)) + ggtitle("Saprotrophs") + labs(x = "Decomposing carcass") + 
  scale_fill_manual(values=c("#00BFC4","#F8766D"))

# SAPROTROPH DIVERSITY # 

tgc <- summarySE(exploration.guild.RA, measurevar="saprotroph.diversity", groupvars=c("decomposing.carcass"), na.rm = TRUE)
tgc

saprotroph.diversity <- ggplot(tgc, aes(x=reorder(decomposing.carcass, -saprotroph.diversity), y= saprotroph.diversity, fill=decomposing.carcass)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() +
  theme_classic() + ylab("Species diversity") + geom_errorbar(aes(ymin=saprotroph.diversity-se, ymax=saprotroph.diversity+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.y = element_text(size = 18), axis.text.y = element_text(size = 18), axis.text.x = element_text(size = 18), axis.title.x = element_text(size = 18),
        legend.position="none", plot.title = element_text(size=20)) + labs(x = "Decomposing carcass") + scale_fill_manual(values=c("#00BFC4","#F8766D")) + ggtitle("Saprotrophs")

# MEDIUM ABUNDANCE # 

tgc <- summarySE(exploration.guild.RA, measurevar="medium.distance", groupvars=c("decomposing.carcass"), na.rm=TRUE)
tgc

tgc$se[2] <- 0.02

medium.distance <- ggplot(tgc, aes(x=reorder(decomposing.carcass, -medium.distance), y= medium.distance, fill=decomposing.carcass)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() +
  theme_classic() + ylab("Relative abundance") + geom_errorbar(aes(ymin=medium.distance-se, ymax=medium.distance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.y = element_text(size = 18), axis.text.y = element_text(size = 18), axis.text.x = element_text(size = 18), axis.title.x = element_text(size = 18),
        legend.position="none", plot.title = element_text(size=20)) + ggtitle("Medium-distance types")+ labs(x = "Decomposing carcass") + 
  scale_fill_manual(values=c("#00BFC4","#F8766D"))

# MEDIUM DIVERSITY # 

tgc <- summarySE(exploration.guild.RA, measurevar="medium.diversity", groupvars=c("decomposing.carcass"), na.rm=TRUE)
tgc

medium.diversity <- ggplot(tgc, aes(x=reorder(decomposing.carcass, -medium.diversity), y= medium.diversity, fill=decomposing.carcass)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() +
  theme_classic() + ylab("Species diversity") + geom_errorbar(aes(ymin=medium.diversity-se, ymax=medium.diversity+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.y = element_text(size = 18), axis.text.y = element_text(size = 18), axis.text.x = element_text(size = 18), axis.title.x = element_text(size = 18),
        legend.position="none", plot.title = element_text(size=20)) + ggtitle("Medium-diversity types")+ labs(x = "Decomposing carcass") + 
  scale_fill_manual(values=c("#00BFC4","#F8766D"))

# long distance ra

tgc <- summarySE(exploration.guild.RA, measurevar="long.distance", groupvars=c("decomposing.carcass"), na.rm = TRUE)
tgc

long.ra <- ggplot(tgc, aes(x=reorder(decomposing.carcass, -long.distance), y= long.distance, fill=decomposing.carcass)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() +
  theme_classic() + ylab("Relative abundance") + geom_errorbar(aes(ymin=long.distance-se, ymax=long.distance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.y = element_text(size = 18), axis.text.y = element_text(size = 18), axis.text.x = element_text(size = 18), axis.title.x = element_text(size = 18),
        legend.position="none", plot.title = element_text(size=20)) + ggtitle("Long-distance types")+ labs(x = "Decomposing carcass")  + 
  scale_fill_manual(values=c("#00BFC4","#F8766D"))

# long distance diversity

# just at Hansen
temp <- exploration.guild.RA[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 11),]

# all data
temp <- exploration.guild.RA[which(exploration.guild.RA$Distance < 11),]
temp <- temp[!is.na(temp$long.diversity),]
tgc <- summarySE(exploration.guild.RA, measurevar="long.diversity", groupvars=c("decomposing.carcass"), na.rm=TRUE)
tgc

long.diversity <- ggplot(tgc, aes(x=reorder(decomposing.carcass, -long.diversity), y= long.diversity, fill=decomposing.carcass)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() +
  theme_classic() + ylab("Relative abundance") + geom_errorbar(aes(ymin=long.diversity-se, ymax=long.diversity+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.y = element_text(size = 18), axis.text.y = element_text(size = 18), axis.text.x = element_text(size = 18), axis.title.x = element_text(size = 18),
        legend.position="none", plot.title = element_text(size=20)) + ggtitle("Long-distance types")+ labs(x = "Decomposing carcass")  + 
  scale_fill_manual(values=c("#00BFC4","#F8766D"))

p <- grid.arrange(saprotroph.richness, saprotroph.diversity, long.ra, medium.distance, medium.diversity, nrow=2, ncol=3)
p

ggsave(plot=p, width=13, height=10, dpi=300, filename = "decomposing.carcass.jpg")



